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ABSTRACT 



We present the results from 9 years of optically monitoring the gravitationally lensed zqso = 0.658 quasar RX Jl 13 1 — 123 1 . The R 
band light curves of the 4 individual images of the quasar are obtained using deconvolution photometry, for a total of 707 epochs. 
Several sharp quasar variability features strongly constrain the time delays between the quasar images. Using three different numerical 
techniques, we measure these delays for all possible pairs of quasar images, while always processing the 4 light curves simultaneously. 
For all three methods, the delays between the 3 close images A, B and C are compatible with being 0, while we measure the delay of 
image D to be 91 days, with a fractional uncertainty of 1.5% (lcr), including systematic errors. Our analysis of random and systematic 
errors accounts in a realistic way for the observed quasar variability, fluctuating microlensing magnification over a broad range of 
temporal scales, noise properties, and seasonal gaps. Finally, we find that our time delay measurement methods yield compatible 
results when applied to subsets of the data. 

Key words. Gravitational lensing: strong - cosmological parameters 
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1. Introduction 

Using the time delays between strong lensing images to measure 
cosmological distances (Refsdal 1964) has a number of advan- 
tages: there is no need for any primary or secondary calibra- 
tor, and there are no effects from the intergalactic or interstel- 
lar medium. The method, originally proposed for gravitationally 
lensed supemovae, has so far exclusively been applied to quasars 
lensed in most cases by individual massive galaxies. Exceptions 
are SDSS J 1004+41 12 and J 1029 +262 3, two quasars lensed by 
galaxy clusters, with long time delays (Fohlmeister et al. 2008 
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. The quasar lens time delay method is now recognized as 
that complements other cosmological probes, in particular 
to constrain Ho as well as the dark energy equation of state pa- 
rameter, w (e.g., |Suyu et al.|2012b[|Linder|201 1 | |Moustakas et al.| 



2009 ). In spite of its advantages, the method has long faced two 



severe limitations to its effectiveness in constraining cosmology. 



* Based on observations made with the 1.2-m Swiss Euler telescope 
(La Silla, Chile), the 1.3-m SMARTS telescope (Las Campanas, Chile), 
and the 1.2-m Mercator Telescope. Mercator is operated on the island 
of La Palma by the Flemish Community, at the Spanish Observatorio 
del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. 
** Light curves will be available at the CDS via anonymous ftp to 
cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz- 
bin/qcat?J/A+A/??? , and on http://www.cosmograil.org. 



First, time delays between the gravitationally lensed images 
of a quasar are hard to measure. The first time delay measure- 
ments were quite controversial (|Vanderriest et aT]|1989 
[eFaLl[T992l |Schild & Thomson|| 19971 IKundic et al.[ 



Understandably, early light curves tended to be short and sparse, 
often too short to clearly demonstrate that microlensing vari- 
ability was not interfering with their analysis. Microlensing is 
seen as an uncorrelated extrinsic variability in the quasar images, 
which results from the time-variable magnification created by 
stars in the lensing gala xy (e.g., |Chang & Refsdal|1979 ; Schmidt 
& Wambsganss 2010]). In the best cases, light curves spanned 
a few years (see, e.g., |Wyrzy kowski et al. 2003} Hjorth et al. 



2002; Burud et al. 2002a]bJ). One consequence is that the numer- 
ical methods used to measure time delays from these light curves 
were exceedingly "optimistic" in their assumptions about extrin- 
sic variability. Later measurements with better data frequently 
yielded delays inconsistent with the error estimates of the earlier 
measurements. 

Second, given measured time delays, and lens and quasar 
image astrometry, there is a famed degeneracy between the time 
delay distance, which is a scale parameter inversely proportional 
to Ho, and the spatial distribution of mass responsible for the 
strong lensing phenomenon. The delays constrain only a com- 
bination of this time delay distance and the surface density of 
the lens near the images ( |Kochan ek 2002). This can be over- 
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come given independent constraints on the structure of the lens. 
Suyu et al. ( |2009 2010) > convincingly show that it is possi- 



ble to control the effects of model degeneracies in the case of 
B 1608+656, a quadruply imaged quasar with accurate radio time 
delays (Fassnacht et al. 2002[ ). To do this, the authors combine 
(1) detailed HST images of the lensed quasar host galaxy, (2) a 
velocity-dispersion measurement of the lens galaxy, and (3) in- 
formation about the contribution of intervening galaxies along 
the line of sight, from galaxy number counts calibrated with nu- 
merical simulations. 

In parallel to the advances in lens modeling, the observa- 
tional situation has drastically evolved as well. Two observa- 
tional groups, COSMOGRAIL and |Kochanek et al.| ( |2006] >, have 
been intensely monitoring * 20 lenses for roughly 10 years. 
In 2010, our 2 groups decided to merge their observational ef- 
forts, with the COSMOGRAIL group focusing on the analysis 
of time delays and the |Kochanek et al.| (j2006 ) group focusing 
on the analysis of microlensing. While preliminary results have 
been published both before and after this merger (Kochanek 
eTaLl[2006] |Vuissoz et aL|2007l [2008) |Morgan et al.|2008a|bl 



exquisite data spanning almost a decade of continuous observa- 
tion are now being released, such as for the quadruply imaged 



2011; Blackburne etaLl 



quasa r HE 0435-1223 ( [Courbin et al 
|20TT) . 

In this paper, we present 9 years of optical monitoring of the 
quadruply imaged quasar RX Jl 131-1231 ( Sluse et al.||2 003 



and measure its time delays with the techniques of |Tewes et a 
(|2012|l. RX Jl 131-1231 is one of the most spectacular lenses of 



our sample. The redshift of the lensing galaxy is zi ens = 0.295, 
while the quasar is at zqso = 0.658. This low quasar redshift 
means (1) that the photometric variations are fast, numerous 
and strong because it is a lower luminosity quasar (see, e.g., 
Mac Leod et al. 12010] ), and (2) that the host galaxy of the lensed 
quasar is seen as a full Einstein ring with many spatially re- 
solved structures in HST images. Similarly, the lensing galaxy 
is sufficiently bright to allow a precise measurement of its ve- 
locity dispersion and possibly of its velocity dispersion profile. 
These characteristics ease both the time delay measurement and 
the lens modeling. The latter, with state-of-the-art inferences of 
cosmological constraints based on our time delay measurements 



of RX J 1 1 3 1 - 123 1 , are presented in |Suyu et aL| ( |2012a| >. 

The COSMOGRAIL and SMARTS observations of 
RX Jl 131-1231 and their reduction are described in Sections [2] 
& [3] while the light curves are presented in Section |4] In 
Section [5] we apply three different curve shifting techniques 
to the light curves and we infer our best measurements of the 
delays along with realistic random and systematic error bars. 
Our results are summarized in Section [6] 



2. Observations 

We have been monitoring the quadruply lensed quasar 
RX Jl 131-1231 (J2000: ll h 31 m 52 s , -12°31'59") since 
December 2003 with 3 different telescopes in the R band (~ 600 
- 720 nm). Table [T] summarizes the observational strategy and 
instrumental characteristics. The light curves presented in this 
paper cover 9 observing seasons (2004 - 2012, see Fig.|4|, with 
707 monitoring epochs in total. The average sampling within the 
seasons is 2.9 days, and the median sampling is 2.0 days. The 
mean seasonal gap in the combined light curves is 132 days with 
a standard deviation of 2 weeks. 

The majority of the measurements for this southern target 
came from the Swiss 1.2-m Euler telescope and the SMARTS 




n 



Euler C2 
Euler ECAM 
SMARTS 1.3-m 
Mercator 




1.0 1.5 2.0 
FWHM [arcsec] 



1.05 1.10 1.15 1.20 1.25 1.30 

e —a/b 



Fig. 1. Distributions of stellar FWHM and elongation e of all 
images contributing to the light curves, by telescope. 



1.3-m telescope, both located in Chile under equally good at 
mospheric conditions. Fig. [T] shows the distributions of the stel 
lar Full Width at Half Maximum (FWHM) and the elongation e 
of each observing epoch, as measured by SExtractor (Bertin & 
|Arnouts|20"T0| > using field stars. The SMARTS 1.3-m telescope 
is guided, in contrast to Euler and Mercator, which are solely 
tracking. This accounts in part for the broad elongation distribu- 
tion of the Euler and Mercator images. 

The original imaging instrument C2 of the Euler telescope 
was replaced in September 2010 by EulerCAM, a liquid nitro- 
gen cooled 4k x 4k e2v 231-84 CCD yielding a pixel scale of 
0'.'215. At the same time the focusing procedure was improved. 
Fig. [T] shows that after two years of observations, images from 
EulerCAM tend to be statistically sharper than C2 images, and 
more comparable to SMARTS 1.3-m data in terms of FWHM. 
All images from the SMARTS 1.3-m telescope were obt ained 
through the optical channel of the ANDICAMM camera ( jDePoy 
|et al.|[2003] i. See |Kochanek eTal] (f2006) for details about the 
SMARTS data. 

Fig. [2] shows the central part of a deep combination of Euler 
monitoring images, totaling 2.5 days of exposure. In the best in- 
dividual exposures the 4 quasar images are well resolved; images 
A and B have the smallest separation of Y!2. 



3. Image reduction 

In this section we describe the reduction procedure leading to the 
light curves of the multiple quasar images of RX J 1 1 3 1 — 123 1 . 
This same procedure will homogeneously be applied to other 
lenses monitored by COSMOGRAIL and SMARTS, hence we 
provide a general description of our methods. 

All exposures are corrected for bias/readout effects and flat 
fielded using standard methods. This prereduction is done in- 
dividually for each telescope and instrument. For Euler and 
Mercator we have developed custom python pipelines that en- 
able efficient human inspection and validation of each step. 
Particular attention is given to a semi-automated selection and 
combination of sky-flats, which is controlled through a simple 
graphical interface highlighting the temporal evolution of instru- 
mental contaminations such as dust. Using difference images we 
check that all sources in the flatfield exposures are adequately 
masked and that the contamination does not evolve significantly 
within sets of flatfields that are stacked as master-flats. The pe- 
riod over which we stack flatfields can be from single days to 
several weeks. Some exposures of the Mercator telescope are 
prereduced using "superflats", combinations of masked science 



http : //www . astronomy . ohio- state . edu/ANDICAM/ 
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Table 1. Overview of our optical monitoring of RX Jl 131-1231 



Telescope 


Location 


Instrument 


Pixel scale 


Field 


Exp. time 


Epochs 


Sampling 


Euler 1 .2 m 


ESO La Silla, Chile 


C2 


0'.'344 


ii' x ir 


5 x 360s 


265 


5.0 (4.0) 


Euler 1 .2 m 


ESO La Silla, Chile 


EulerCAM 


0'.'215 


14' x 14' 


5 x 360s 


76 


5.4 (4.1) 


SMARTS 1.3 m 


CTIO, Chile 


ANDICAM 


0'.'369 


6.1' x 6.1' 


3 x 300s 


288 


6.1 (5.1) 


Mercator 1 .2 m 


La Palma, Canary Islands, Spain 


MEROPE 


0'.'193 


6.5' x 6.5' 


5 x 360s 


78 


8.9 (4.5) 



Notes. The column "Exp. time" indicates the number of dithered exposures per epoch and their individual exposure times. One epoch corresponds 
to one data point in each light curve. The temporal sampling of the observations is given as the average (median) number of days between 
consecutive epochs, excluding the seasonal gaps. 
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Fig. 2. Part of the field of view of the Swiss 1.2-m Euler telescope around RX Jl 1 3 1 — 123 1 . The wide-field image is a combination 
of 600 exposures of 360s each, corresponding to a total exposure time of 2.5 days. The inset shows a single 360s exposure from the 
EulerCAM instrument under good conditions (FWHM~1", e = 0.1, no Moon). All images are taken in the R band. We mainly use 
the stars labeled 1 to 4 to build a PSF model for each exposure. Stars 1 to 5 are used for the photometric calibration. 



exposures obtained in a single night. We have developed vari- 
ants of these pipelines adapted to the particularities of other 
COSMOGRAIL telescopes. 

From here on, all images are processed by a single deconvo- 
lution photometry software package. The principle ideas for this 
reduction procedure are the same as employed for previously 
published COSMOGRAI L light curves ( |Vuissoz et aE1[2007 



I Courbin et al. 2011 



). Over the years we have reworked the 
steps, and implemented the procedure in the form of a python 
pipeline linked to a relational database, containing one entry per 
exposure. 



3.1. PSF construction 

We first build a conservatively smooth sky model, obtained using 
SExtractor, for each exposure. This sky model is not critical, as 
the subsequent photometry procedure will fit for any residual sky 
level around all sources of interest. 

We then align the images, separately for each instrument, 
and individually estimate the Point Spread Function (PSF) of 
each exposure. This is done by fitting several field stars using a 



common model, composed of (i) a simply parametrized profile 
and (ii) a regularized fine pixel array. The details of this proce- 
dure, which is part of a general purpose deconvolution package, 
are described in Cantale et al. (in prep). For most of the expo- 
sures of RX J1131-1231,we use the stars labeled 1 - 4 in Fig. [2] 
to build the PSF model. The pipeline allows us to easily explore 
and compare different choices of PSF-stars. The final choice of 
stars is empirically selected so as to yield the least scatter in the 
final quasar light curves. In the case of RX Jl 131-1231 the sit- 
uation is close to optimal, since (1) the stars 1 - 4 are bright 
but still in the linear regime of the CCDs, (2) they surround the 
lens at modest angular distances, and (3) all stellar companions 
or background objects can be identified and masked, yielding 
clean fitting residuals. Cosmic ray hits and CCD artifacts are au- 
tomatically masked using a variant of the L. A. Cosmic algorithm 
(van Dokkum 2001 ). We visually supervise the PSF construction 



and manually adapt the selection of PSF-stars for problematic 
frames. The construction of the PSF models is the predominant 
technical issue affecting the final quality of the light curves. 
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3.2. Photometric normalization of the exposures 

We next compute a multiplicative normalization coefficient for 
each exposure, based on photometry of field stars. These coeffi- 
cients will constitute the reference for the differential photome- 
try of the quasar images. 

In a first step, this computation is done separately for each in- 
strument. We measure the instrumental fluxes (in photons) 
of star i in each exposure j, by fitting the PSF ; of this exposure 
to each star, leaving only the fluxes, astrometric shifts, and back- 
ground levels of each star as free parameters. For each star and 
exposure, we then compute an individual calibration coefficient 



|et al.|2003), as well as in past COSMOGRAIL publications (e.g., 



med/A^y) 



(1) 



where medj denotes the median over all exposures of this instru- 
ment. Then, for each exposure j, we obtain the normalization 
coefficient through 



med,(c;p. 



(2) 



The advantages of this simple and fast method is that it does 
not select a single exposure as the "reference" for all stars. For 
each instrument, the distribution of these coefficients Cj is highly 
unimodal, as the exposure time is kept constant. 

To select the normalization stars we first inspect the normal- 
ized stellar light curves for any anomalies or variability on scales 
of months or years, and, if required, we iteratively repeat the co- 
efficient computation with a refined selection of stars. In terms of 
the smoothness of the final quasar light curves, the best results 
are generally obtained by computing these coefficients from a 
few well exposed stars with a comparable and high signal-to- 
noise ratio, as opposed to using larger numbers of noisier stars. 
Furthermore, the ideal normalization stars should be of similar 
color to the quasar. This minimizes potential differential effects 
due to varying atmospheric conditions given the relatively broad 
R band filters used for the monitoring. 

In a second step, we rescale the coefficients of each instru- 
ment so that a star whose color is closest to the quasar gets the 
same normalized median flux across all instruments. Residual 
adjustments to this normalization between instruments will be 
performed later, using the quasar light curves in temporal re- 
gions where the data from different telescopes overlap. 

3.3. Photometry of the quasar images 

We obtain deblended light curves of the quasar images by 

1998). 



MCS deconvolution photometry, following Magain et al. 



This algorithm fits a single model consisting of (1) some num- 
ber of point sources and (2) a fine pixel array (hereafter the 
pixel channel) simultaneously to all exposures. In the case of 
RX Jl 13 1 — 1231, 4 point sources model the quasar images, while 
the pixel channel represents both the lensing galaxy and the 
lensed host galaxy (Einstein ring). This model is scaled by the 
normalization coefficients from Equation [2] and convolved by 
the exposure-specific PSFs, before being fit to the data. By con- 
struction, the relative astrometry of the 4 point sources and the 
pixel channel are common to all exposures. Only the fluxes of 
the point sources, the absolute astrometric shift, and a flat resid- 
ual sky level around the lens are free to vary between expo- 
sures. Iteratively, all these parameters are optimized, together 
with the structure of the regularized pixel channel, so as to min- 
imize a single global x 1 ■ This algorithm was already success- 
fully applied in the discovery paper of RX Jl 131 — 1231 (Sluse 



Eigenbrod et al.|2007[ |Vuissoz et al.|2008| |Courbin etal.|2011 
Slus e et al.||2012|> or for s imilar monitoring data (Burud et al. 
2002a||Hjorfh et al.|2002}|Jakobsson et al.|2005||Morgan etaL 
|2012| ). The light curves of the quasar images are a direct output 
of this procedure. 

Errors in the astrometry or the structure of the pixel channel 
might degrade or bias the photometry of the quasar images. We 
observe that the influence of the astrometry on the light curves 
is small. Even for bright quasar images, alterations of the rel- 
ative position of the point sources by as much as 0'.'05 do not 
significantly modify the light curves. Larger errors tend first to 
smoothly bias the magnitude measurements, before introducing 
additional scatter. To obtain light curves for gravitational lenses 
with image separations on the order of the resolution of the best 
monitoring data, we find no difference between using the tight 



astrometric constraints from HST images (Mor gan et al.[ :2006 



Chantry et al. 2010) or letting the deconvolution algorithm freely 
optimize the astrometry of the model. 

We reach similar conclusions about the impact of the pixel 
channel, and its mandatory regularization. Simple numerical ex- 
periments show that to first order, the effect of different plausible 
solutions is similar to very small additive shifts to the flux of the 
quasar light curves, with only marginally perceptible effects even 
for faint quasar images. In practice we systematically constrain 
this pixel channel as well as the point source astrometry using 
only a subset of images with the best resolution. 

3.4. Photometric error estimation 

In this section we describe how we compute a rather formal "best 
case" error estimate for the photometry of each quasar image in 
each exposure. The normalized flux of a quasar image in a given 
exposure can be written /* = A* • c, where A* is the measured 
number of photons, and c is the normalization coefficient from 
Equation [2] We assume that the random error on /* has two in- 
dependent sources: (i) the shot noise cr Nt of the quasar image 
itself, and (ii) the noise cr c of the normalization coefficient as 



computed in Section 3.2 



We compute cr N ^ following the standard "CCD equation" 
(see e.g. |Howell|2006l chap. 4.4) 



+ » pix • (S + R 2 ) 



(3) 



where S is the sky level, R is the CCD read noise (both in pho- 
tons per pixel), and n p ; x is the number of (equivalent) pixels of 
the software aperture. The dark current is negligible in our im- 
ages. MCS deconvolution of point sources corresponds to PSF 
fitting, and following Heyer & Biretta (2004 section 6.5.1) we 
use for n pix the reciprocal of the PSF sharpness, 



— = V (PSF,,,,) 2 , 



(4) 



where PSF/,„ is the fraction of light in the total PSF at pixel Im. In 
this simple "best case" computation, we assume that the PSF is 
perfectly known, and that the sky level is well determined by the 
fitting procedure. We do not take into account the surface bright- 
ness of the lens and host galaxies, which is generally negligible 
with respect to the sky level. Furthermore we do not compute 
noise terms due to the blending with other point sources. 

As a sanity check, we have tested our deconvolution photom- 
etry algorithm on simulated point source images created with 
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SkyMaker ( |Bertin|2009| ) for various realistic signal-to-noise ra- 
tios, elongations, and seeing conditions. Using two or more well- 
exposed stars to build the PSF, the statistical scatter obtained by 
our photometry pipeline is only marginally larger than the best- 
case shot noise computed through the above procedure. Using 
similar synthetic images, we have also verified that deconvolv- 
ing two partially blended sources of different fluxes does not sig- 
nificantly bias their measured flux difference in any direction. 

To estimate the random uncertainty of the normalization co- 
efficient, we use the standard error of the mean (SEM) of the 
individual coefficients obtained from the n normalization stars 
in the exposure (see Equation[T} 



<t c = 



where s' 



(5) 



Finally we obtain the error estimated ct^ccd on the normalized 
flux through 



N: 



(6) 



u /„ CCD 

IT 

Note that for our COSMOGRAIL and SMARTS monitoring 
data, the shot noise term from the CCD equation dominates this 
error budget. The contribution from o~ c is typically <k 3% of 
o"/, ccd- Such small calibration errors are important only when 
measuring very short delays, where they would introduce posi- 
tively correlated noise into the curves. 

3.5. Combination of points per night 

In each monitoring night, we observe the lenses m times (m — 3 
or 5) over * 30 minutes, yielding flux measurements fa, with 
k — 1, .. , m for each quasar image. To reduce the CPU cost and 
to reject outliers, we bin these measurements by epochs, sepa- 
rately for each instrument and telescope, before measuring the 
time delays. We attribute to each epoch the medians of the pho- 
tometric measurements fa within that night, and the mean of 
the Heliocentric Julian Dates (HJD). 

This approach also allows us to obtain an empirical photo- 
metric error estimate for each quasar image, using a measure of 
the spread of the image's fa- For increased robustness against 
outliers, we quantify this spread using the Median Absolute 
Deviation from the median, or MAD (Hoaglin e t al.|1983 1, 

MAD(A*) = med( | fa - med(fa) | ). (7) 

To estimate the usual cr of an (assumed) Gaussian distribution, 
the MAD is rescaled as 



07. mad = 1-4826 -MAD(fa). 



(8) 



These procedures give us two different error estimates for each 
epoch and for each quasar image: (1) the median of the errors es- 
timated individually for each exposure med(o"y t ccdX an d (2) the 
more empirical cr^ MAD . We observe, as expected, that both er- 
ror estimations are highly correlated, but also that the empirical 
07, mad is typically twice as big. 

The uncertainty we finally assign to each epoch and quasar 
image uses the larger of med(cry tCCD ) and tr^ mad- Indeed a real- 
istic error estimate cannot be smaller than either of them. We di- 
vide this estimate by the square root of the number of exposures 
m of the epoch, so that the final error estimate for the median 
photometric measurement on a given instrument and night is 



r med(/ t ) 



max[med(o- /tCCD ), 07; mad] 

-0M 



Euler C2 
Euler ECAM 
SMARTS 1.3-m 
Mercator 




19 18 
R [mag] 

Fig. 3. Photometric lcr error estimates of each observing epoch 
as a function of approximative R band magnitude. These errors 
are estimated following Equation [9] 

Fig.[3]displays the distribution of the resulting error estimates as 
a function of R band magnitude and instrument. The data points 
from all 4 quasar images of RX Jl 131-1231 are shown, yielding 
a broad spread in magnitudes. 

3.6. Combination of telescopes 

The lens galaxy and the lensed images of the quasar host galaxy 
differ in color from the quasar itself. As a consequence, the dif- 
ferent R filters and CCD response functions of the monitoring in- 
struments might "see" these contaminating sources with slightly 
different amplitudes even with a perfect calibration of the quasar 
fluxes. This can result in small mismatches between the light 
curves obtained using different combinations of telescopes and 
instruments. 

We correct our light curves for such small effects (occasion- 
ally as large as 10% of the flux of the faintest quasar images) in 
this final step. We typically select the instrument with the longest 
or the highest quality curve as a reference. For RX Jl 131-1231, 
this is the SMARTS light curve. For each of the other instru- 
ments, we optimize additive magnitude and flux shifts (i.e., mul- 
tiplicative and additive flux corrections) so to minimize a dis- 
persion measure between each instrument's light curve and the 
reference light curve. We compute this dispersion following the 
curve shifting technique presented in |Tewes et al.| {2012), but 
evaluate it between the light curves of different instruments, in- 
stead of different quasar images. Provided that the colors of the 
quasar images are not differentially reddened by absorption in 
the lens galaxy, we optimize a single common magnitude shift 
per instrument, and individual flux shifts for each quasar image 
and instrument. This is adequate for RX Jl 13 1 — 123 1, as previ- 
ously suggested by Sluse et al. (2007} and Chartas et al. ( |2009 
|20L2l >. 



4. The light curves of RX J1 131 -1231 

All these steps were applied to the COSMOGRAIL and 
SMARTS data for RX J 1 1 3 1 - 123 1 . Fig.g]shows the final 9-year 
long light curves of the 4 quasar images. This data is available in 
machine-readable form at the CDS and on the COSMOGRAIL 
websit^] The light curves are dominated by intrinsic quasar 
variability, with some features on scales as fast as a few weeks. 
It can be readily seen in the 2008 season for instance, that the 
delays between A, B and C must be very small, while D is de- 
layed by a bit less than 100 days. Intriguingly, looking only at 



(9) 



http : //www . cosmograil . org 
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Fig. 4. Optical monitoring of RX Jl 13 1 — 123 1 , as obtained from deconvolution photometry. From top to bottom are shown the 
R band light curves for the quasar images A, B, C, D along with the Icr photometric error bars. Colors encode the contributing 
instruments. Curves B and D have been shifted by +0.3 magnitudes for display purposes. The light curves are available in tabular 
form from the CDS and the COSMOGRAIL website. 



the first season of A, B and C, one might guess that A is signif- 
icantly delayed with respect to B and C. We attribute this dis- 
crepancy to microlensing variability, which manifestly changes 
the magnitude difference between the A and B images from the 
first to the second season. We discuss this "event" further in 
Section |5.3| Prominent microlensing variability on large time 
scales is omnipresent, as the flux ratios between the quasar im- 
ages evolve by as much as a magnitude. These microlensing ef- 
fects in RX Jl 13 1 — 123 1 have been analyzed in [Morgan et al.| 
< |20T0] ) and |Daietal.| ( |20T0| . 

Lastly, we observe that the photometric error estimates, ob- 
tained from equation [9] match the observed scatter well in the 
smooth sections of curves from the individual telescopes. They 
are certainly not conservatively large, but we stress that the scale 
of these error estimates has no direct influence on the uncertain- 
ties that we will compute for the time delay measurements in the 
next Section. Our results are robust against a deliberate increase 
of these error estimates by up to a factor of 5. 



5. Time delay estimation 

In this section we infer the time delays of RX J 1 13 1 — 123 1 
from the light curves shown in Fig. [4j closely following the 
curve shifting and uncertainty evaluation procedures described 
in Tewes et al. (2012 1. We summarize the principal ideas be- 
low. A major difficulty, and potential source of bias, for curve 
shifting methods is the presence of "extrinsic" variability in the 
light curves, on top of the "intrinsic" quasar variability common 
to all 4 images. The main source of the extrinsic signal is vari- 
able microlensing magnification due to the motions of the stars 
in the lens galaxy. As shown by |Mosquera & Ko chanek (2 01 l| l, 
microlensing can affect light curves over a broad range of time 
scales. For RX J 1 1 3 1 - 123 1 , [Mosquera & Kochanek| ( pOTT) es- 
timate a time scale of » 1 1 years for the crossing of a stellar 
Einstein radius, and « 3 months for the source radius to cross a 
caustic. 

The curve shifting methods of |Tewes et al.| ( [20T2] ) try to min- 
imize the bias due to such extrinsic variability in different ways. 
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They all rely on iterative optimizations of time shifts of the 4 
light curves, and yield self-consistent point estimates of the time 
delays between all 6 image pairs : 

1. The free knot spline technique simultaneously fits one 
common intrinsic spline for the quasar variability, and in- 
dependent, smoother extrinsic splines for the microlensing 
to the light curves. The curves are shifted in time to optimize 
this fit. This method is similar to the polynomial method of 
|Kochanek et"aT1 ( |2006l >. 

2. The regression difference technique shifts continuous re- 
gressions of the curves in time, to minimize the variability 
of the differences between them. This novel method does not 
involve an explicit model for the microlensing variability. 

3. The dispersion-like technique, inspired by |Pelt et aT 
( |1996] ), shifts the curves so as to minimize a measure of the 
dispersion between the overlapping data points. This method 
has no explicit model for the common intrinsic variability of 
the quasar, but it includes polynomial models for the extrin- 
sic variability. 

Using several independent algorithms allows us to cross- 
check our estimates for technique-dependent biases. Indeed, as 
important as the time delay point estimates themselves is the 
reliable estimation of their uncertainties. For this we follow a 
Monte Carlo approach, by applying the curve shifting techniques 
to a large number of synthetic light curves with known time de- 
lays. These curves are drawn from a light curve model that mim- 
ics both the observed intrinsic and extrinsic variability of the 
real observations, randomizing only the unrecoverable short time 
scale extrinsic variability. This latter "correlated noise" locally 
adapts its amplitude to the scatter of the observed data points. 
As a result, the synthetic curves are virtually indistinguishable 
from the real ones (see Figs. 5 and 6 of Te wes et al.|2012} for an 
illustration). 
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Fig. 5. Analysis of the "intrinsic variance" of the time delay 
point estimators, found by running the three curve shifting tech- 
niques 200 times on the light curves shown in Fig. [4] starting 
the optimizations from random initial time shifts. The widths of 
these distributions reflect the failure of the methods to converge 
to a single optimal solution given the data. We stress that these 
distributions do not represent probability density functions the 
time delays, and that a small intrinsic variance does not imply 
that a time delay estimation is precise or accurate, only that the 
error surface is relatively smooth. 



5.1. Applica tionto RXJ1 131-1231 

We begin by evaluating the robustness of the time shift optimiza- 
tion of each method given the data. For this we repeatedly run 
the delay point estimators on the observed light curves. Each 
time, we start the optimizations from random initial conditions, 
using time shifts uniformly drawn ±10 days around plausible so- 
lutions. The resulting distributions of point estimates are shown 
in Fig. |5J characterizing what we call the intrinsic variance of 
the estimates. 

These tests are essential to check the reliability of the non- 
linear optimization algorithms for a given data set and model pa- 
rameters. We stress that this procedure should not be interpreted 
as importance sampling of any posteriors. Generally speaking, 
overly flexible models dilute the time delay information and 
yield higher intrinsic variances. For the free knot spline tech- 
nique we have chosen an average knot step of 20 days for the 
quasar variability spline, and 150 days for the extrinsic splines. 
For the dispersion-like method we model the extrinsic variability 
by independent linear trends on each season. If the light curves 
sufficiently constrain the time delays, the free knot spline and re- 
gression difference techniques can easily be adjusted to display 
small intrinsic variance and roughly unimodal distributions, as 
in Fig. [5] The dispersion-like technique is inherently more sen- 
sitive to initial conditions, due to a much higher roughness of 
the scalar objective function that it minimizes. As discussed in 



Tewes et al. ( 2012 1, it is important to note that smoothing the ob- 



jective function does not necessarily lead to more accurate time 
delay estimates. For each technique, we use the mean values of 



the distributions in Fig.[5]as our best time delay estimates. They 
correspond to the points in Fig. [6] which summarizes our results 
forRX J1131-1231. 

The remaining part of the analysis is solely about estimating 
realistic error bars for these point estimates. For this we proceed 
by blindly running the exact same curve shifting techniques on 
1000 sets of fully synthetic light curves with true time delays 
randomly chosen around the measured ones. Again, we start the 
methods from random initial shifts. Statistics of the resulting 
time delay measurement errors (i.e., measurement - truth), are 
shown in Fig. [7] as a function of the true delays of the synthetic 
curves. In this figure, the shaded rods show the mean measure- 
ment error, which we call systematic error or bias, while the error 
bars indicate the standard deviation of the measurement errors, 
representing the random errors. The total errors for our delay 
measurements, shown by the error bars in Fig. [6] are computed 
by adding the maximum bias and the maximum random error 
of each panel in quadrature. Note that we run the curve shift- 
ing techniques only once on each random synthetic curve set, 
so this error analysis takes into account the intrinsic variance. 
This is rather conservative, given that we use the mean result 
of several time shift optimizations as our best estimates for the 
observed data. However, the additional contribution to the un- 
certainty estimates is negligible for methods with low intrinsic 
variance. Before discussing these results, we show in Fig.[8]the 
same measurement errors as in Fig. [7] but plotted for each delay 
against each other, in order to explore potential abnormal corre- 
lations. 
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Fig. 6. Time delay measurements along with the lcr errors for RX Jl 13 1 — 123 1 obtained by our 3 standard techniques from the 
full 9-season long light curves shown in Fig |4] The random and systematic error contributions are given in parentheses for each 
delay. The error bars represent the random plus systematic errors summed in quadrature. A positive AB-delay A?ab means that 
image B leads image A. We consider the measurements from the regression difference technique, which display the smallest bias 
and variance in our error analysis, as the delays that should be used to constrain cosmology and lens models. 
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Fig. 7. Results of the Monte Carlo analysis leading to the error estimates for our time delay measurements shown in Fig.|6| We 
obtain our uncertainty estimates by applying the curve shifting techniques to 1000 synthetic light curve sets that closely mimic the 
observed data but have known "true" time delays. The vertical axes show the delay measurement error, as compared to the true 
delays used to generate the synthetic curves (horizontal axes). Separately for each panel, the outcomes are binned according to 
the true time delays. The bin intervals are shown as light vertical lines. Within each bin, the shaded rods and error bars show the 
systematic and random errors, respectively, of the delay measurements for each technique. 
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Fig. 8. Correlations of delay measurement errors for synthetic 
light curves mimicking RX J 1 1 3 1 — 123 1 . The measurement er- 
rors are the same as shown in Fig. [7J but this time marginaliz- 
ing over the true delays. Crosshairs indicate zero error, and the 
inset shows the scale of each panel. For clarity, only single con- 
tours at half of the maximum density are shown for two of the 
techniques. We observe no correlations (i.e., oblique contours) 
between the unrelated delay measurements along the short diag- 
onal of this figure. 



5.2. Discussion 

In terms of simple lens model considerations, RX Jl 131-1231 
is a "long axis cusp lens", with the source is located inside a cusp 
of the tangential caustic curve on the long axis of the potential 
(Slus e et al.|2003| >. Images A and D are the saddle points of the 
arrival time surface, while B and C are minima ( |Blandford &| 
Narayan|T9 86). Fig. 14 of Sa ha et al.| (|2~006 ) gives an illustration 
of this particular lens. Hence the delays A?ab and At^c are pre- 
dicted to be small but positive, while Atpj) is negative and large; 
the possible arrival time orders are BCAD or CBAD. 

Our measured delays, as shown in Fig. [6] are consistent with 
these predictions, although the delays between A, B and C are 
compatible with being zero given the lcr errors. It is reassuring 
to observe that the three methods yield consistent results for all 
delays despite the very different methodologies, which we see as 
a success of our error estimation procedure. Keep in mind how- 
ever that the delay measurements are not independent, because 
they all use the same single set of observed light curves. 

Let us look in more detail at the measurement error analy- 
sis shown in Fig. [7] which is based solely on synthetic curves. 
Overall, we do not observe any strong dependence of a method's 
bias or random error on the true time delays used in the simula- 
tions. For most image pairs and methods, the random errors are 
more important than the almost negligible bias. Exceptions are 
the two results from the dispersion-like technique, which is ob- 
served to overpredict At^B, and underpredict Atsc, by about one 
day. Remarkably, this same technique also measured the high- 
est (lowest) delay At^B (A?bc) f° r the observed curves (Fig. [6j, 
when compared to the other methods. We see this as a further in- 



dication that the bias estimates are reliable. Finally, Fig.|8]shows 
no sign of abnormal correlations between measurement errors of 
quasar image pairs, whether they share a common quasar image 
or not. By construction our time delay uncertainty estimates for 
each image pair marginalize over all other pairs. 

Which delay measurement method performs best on 
RX J 1 1 3 1 — 123 1 ? Given that the regression difference technique 
yields both the smallest biases and the smallest random errors, 
we simply conclude that its measurements are the most informa- 
tive ones. We will use the delays from the regression difference 
technique, expressed with respect to quasar image B, to mea- 
sure the time-delay distance towards RX J 1 13 1 — 123 1 . Details 
and results of lens modeling, as applied to RX Jl 131-1231, are 
described in detail in Suyu et al. ( 2012a| l. 



5.3. Delay measurement on subsets of the light curves 

Our long light curves, featuring several delay-constraining in- 
trinsic variability patterns, suggest that we should independently 
measure the time delays from subsets of the observing seasons. 
This analysis represents an invaluable check of the consistency 
of the time delay estimation procedure of |Tewes et al. ( 2012| l, 
and hence for the results from the full light curves. 

We analyze 3 subsets of the available data : (1) the first two 
seasons, (2) the first 5 seasons, and (3) the remaining last 4 
seasons. Case (1) is chosen for comparison with |Morgan et al7| 
(2006). We perform the analyses from scratch, without using 
any knowledge about the delay estimates from the full curves. In 
particular, for each case, we build a new model for the synthetic 
light curves, independently adjusted so that they best resemble 
the observations in that time period based on the statistical cri- 

pOTIl Section 7). We do not alter 



teria presented in Tewes et al. 



any parameters of the curve-shifting techniques. 

Fig. [9] presents the resulting delay measurements. The data 
points and error bars depict the individual point estimates and 
the corresponding lcr total errors obtained by each technique for 
the 3 cases. The shaded regions show the lcr intervals from the 
full 9 seasons taken from Fig. [6] We observe : 



Using only the first 2 seasons, our three methods are system- 
atically biased towards large values of At^B and At^c- Our 
interpretation is that some microlensing variability in image 
A conspiratorially imitates a time shift, particularly around 
mid 2004. The estimates from the dispersion-like and spline 
technique roughly reproduce the results obtained from the 
same two seasons by Morgan et al. (2006) At^B = 12.0+ 1.5, 
At AB = 9.6 ±1.9 and AfAzj = -87 ± 8, but with significantly 
wider, yet still too small, error bars. This demonstrates the 
need of long monitoring programs to measure accurate time 
delays and minimize the influence of unfortunately placed 
microlensing events. 

For these same two seasons, the regression difference 
method yields far better estimates for these delays - includ- 
ing adequately sized error bars that encompass the delays 
estimated from the full light curves. Note that the regression 
difference technique is the only one that does not involve an 
explicit model for the microlensing. 

For the other 2 divisions of the data, namely seasons 1-5 and 
6-9, which are disjoint and hence independent, our methods 
yield consistent results. As would be expected, both subcases 
show larger error bars than the combined analysis of all sea- 
sons. 
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Fig. 9. Application of the curve shifting techniques to subsections of the full light curves. The square diamonds (top) show mea- 
surements using only the 2 first seasons, the leftward triangles (middle) use the first 5 seasons, and the rightward triangles (bottom) 
use the last 4 seasons. All error bars depict total errors, as in Fig. [6] The blue and green shaded regions correspond to the interval 
covered by the total error bars using the full 9 seasons, by the spline and regression difference techniques, respectively. 



6. Conclusions 

The first part of this paper describes the COSMOGRAIL data 
reduction procedure, which will be used to reduce all the 
data gathered by our monitoring campaign of gravitationally 
lensed quasars. In the second part we apply this pipeline to our 
COSMOGRAIL and SMARTS observations of the quad lens 
RX Jl 131-1231, leading to an unprecedented set of 9 year long 
light curves of high photometric quality. Several strong and fast 
intrinsic quasar variability patterns constrain the time delays 
between the multiple images. Microlensing-related "extrinsic" 
variability is clearly present, as pointed out and analyzed in pre- 
vious studies (|Sluse et al |2006| |2W7t |Morgan et al.|2010[|Dai 



and (2) the large scale structures along the line of sight to the 
quasar ( Bar-Kana|1996 1. 



et al.|2010 Chartas et al.|2012 1. However, this distorting signal 
does not prevent us from measuring accurate time delays, using 
the three independent algorithms of |Tewes et al. ( 2012| l. 

The best time delay estimates of RX J 1 1 3 1 — 123 1 are pro- 
vided by the regression difference technique. It measures the 91- 
day delays between D and the other quasar images to a frac- 
tional lcr uncertainty of 1.5%. This error estimate is obtained 
by applying the techniques to synthetic curves with known time 
delays, which contain extrinsic variability features similar to 
the observed ones. We demonstrate the consistency of our er- 
ror estimates by independently measuring time delays - includ- 
ing error bars - from subsets of the observed light curves of 
RX Jl 131-1231. This experiment also reveals that long multi- 
year monitoring is essential to reliably measure time delays, de- 
spite progress on the methods. 

The results from this paper are used to constrain the time 
delay distance towards RX Jl 131-1231, and deduce stringent 
implications for cosmology in |Suyu et aL] ( |2012a| l. With our 
time delay measurement errors of only 1.5%, the accuracy of 
this cosmological inference is actually limited by the residual 
uncertainty of (1) the gravitational potential of the lens galaxy 
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